home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Info-Mac 4
/
Info_Mac IV CD-ROM (Pacific HiTech Inc.)(August 1994).iso
/
Science
/
RLaB
/
toolbox
/
erf.r
< prev
next >
Wrap
Text File
|
1994-04-25
|
2KB
|
106 lines
//-------------------------------------------------------------------//
//
// Syntax: erf ( A )
// erfc ( A )
// erfcc ( A )
// inverf ( A )
// Description
// Error function, complimentary error functions, and inverse error
// function.
// erf (A) = 2/sqrt(pi) sum exp(-t^2) (from 0 to A)
// erfc (A) = 1 - erf (A)
// erfcc (A) = erfc (A) (but faster)
// See Also: gamma
//-------------------------------------------------------------------//
// Dependencies
rfile gamma
erf = function ( X )
{
local (e, i);
e = zeros (size (X));
for (i in 1:X.n)
{
if (X[i] < 0)
{
e[i] = -gammp (X[i].^2, 0.5);
else
e[i] = gammp (X[i].^2, 0.5);
}
}
return e;
};
erfc = function ( X )
{
local (e, i);
e = zeros (size (X));
for (i in 1:X.n)
{
if (X[i] < 0)
{
e[i] = 1.0 + gammp (X[i].^2, 0.5);
else
e[i] = gammq (X[i].^2, 0.5);
}
}
return e;
};
//
// A simpler, and much faster version of erfc().
//
erfcc = function ( X )
{
local (ans, i, t, z);
z = abs (X);
t = 1 ./ (1 + 0.5*z);
ans = t.*exp (-z.*z-1.26551223+t.*(+1.00002368+t.*(0.37409196+...
t.*(+0.09678418+t.*(-0.18628806+t.*(0.27886807+...
t.*(-1.13520398+t.*(+1.48851587+...
t.*(-0.82215223+t.*0.17087277)))))))));
for (i in 1:X.n)
{
if (X[i] < 0)
{
ans[i] = 2-ans[i];
}
}
return ans;
};
inverf = function ( X )
{
local (b, bp, f, fd, m, n, tol, tp);
//
// Newton-Raphson iteration on ERF
//
tol = 0.00001;
tp = 1/sqrt (2*pi);
m = X.nr;
n = X.nc;
bp = ones (m, n);
b = zeros (m, n);
while (any (any (abs (b-bp) > tol)))
{
bp = b;
f = ((1-erfcc (bp)) - X)/2;
fd = exp (-0.5*bp.*bp).*tp;
b = bp - f./fd;
}
return b;
};